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Abstract 


A  method  for  using  integrating  detectors  to  estimate 
the  phase  of  an  optical  wavefront  is  investigated.  The 
phase  and  amplitude  are  assumed  to  vary  slowly  compared 
to  the  integration  time  and  the  integrated  samples  are 
shown  to  be  corrupted  by  white  gaussian  noise.  A  maxi¬ 
mum  a-posteriori  nonlinear  Kalman  filter  is  derived  and 
simulated  for  both  constant  and  random  walk  phase  pro¬ 
cesses.  The  performance  of  the  filter  is  compared  to  its 
Cramer-Rao  lower  bound  and  a  first  order  linearized  phase- 
locked  loop  (PLL).  The  filter  never  performs  a  great  deal 
better  than  the  PLL,  but  it  can  be  implemented  in  a  low 
bandwidth  system  whereas  the  PLL  assumes  a  wideband  system. 

The  local  oscillator  is  initially  assumed  to  be  sta¬ 
bilized  in  frequency  and  amplitude,  and  later  a  frequency 
shift  is  introduced.  The  filter  manages  to  acquire  and 
track  the  phase  in  a  high  carrier-to-noise  ratio  (CNR)  en¬ 
vironment  although  it  does  not  estimate  the  frequency  shift 
well.  At  30  db  CNR,  the  phase  is  estimated  within  about 
0.02  radians,  whereas  the  frequency  shift  estimation  error 
is  on  the  order  of  500  kHz  for  a  frequency  shift  of  2MHz. 
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OPTICAL  PHASE  ESTIMATION 


FROM  INTEGRATED  SAMPLES  OF  THE 
HETERODYNED  WAVEFRONT 

I .  Introduction 

The  phase  of  an  optical  wavefront  cannot  be  measured 
directly  by  a  detector  since  this  would  require  an  extra¬ 
ordinarily  fast  detector  and  associated  processing  elec¬ 
tronics.  All  presently  used  optical  detectors  respond  only 
to  the  incident  total  optical  power,  not  the  instantaneous 
value  of  the  E-field.  Nevertheless,  optical  phase  can  be 
measured  indirectly  via  several  methods,  most  involving 
some  sort  of  heterodyning  system  which  moves  the  carrier  co¬ 
herently  to  a  lower  frequency  where  either  the  spatial  or 
temporal  characteristics  of  the  resulting  interference  pat¬ 
tern  can  be  examined  (Ref .4 : 3301) .  Among  other  things,  op¬ 
tical  phase  can  be  used  to  determine  refractive  index  vari¬ 
ations,  examine  transparent  objects,  perform  nondestructive 
testing  of  surfaces,  and  reconstruction  of  an  incident  light 
wavefront  <for  instance,  to  correct  aberrations  in  an  active¬ 
ly  compensated  optical  system,  whether  introduced  by  the  op¬ 
tics  or  the  atmosphere). 

As  is  demonstrated  in  the  theoretical  development  later 
in  this  paper,  the  output  of  a  detector  in  an  optical 


receiver  is  at  the  heterodyne  frequency,  typically  on  the 
order  of  tens  of  megahertz  or  more.  In  order  to  pass  this 
waveform,  the  detector  must  have  a  bandwidth  of  at  least 
that  much,  and  usually  more  if  doppler  shifts  may  be  en¬ 
countered.  Such  wideband  detectors  are  expensive  and  not 
readily  available  in  large  array  structures  with  uniform 
response  characteristics.  On  the  other  hand,  large  uniform 
arrays  of  charge  coupled  devices  (CCD’s)  are  readily  avail¬ 
able.  These  devices  do  not  have  sufficient  bandwidth  to 
produce  an  output  which  is  directly  proportional  to  the  het¬ 
erodyned  signal  and  thus  produce  a  vector  of  outputs  that 
are  samples  of  the  time  integral  of  the  heterodyned  wave- 
front.  That  is,  CCD's  are  photon  counting  detectors. 

Interestingly,  Wyant  has  shown  (Ref . 20 : 2624 )  that,  in 
the  absence  of  noise,  it  is  possible  to  determine  the  phase 
of  a  wavefront  from  these  integrated  samples.  The  purpose 
of  this  paper  is  to  extend  and  generalize  Wyant's  analysis 
and  to  analyze  how  well  a  system  performs  in  the  presence  of 
corruption  by  white,  gaussian  noise. 

In  this  paper,  the  phase  of  an  optical  field  is  defined 

as  the  argument  of  the  received  optical  field  phasor, 
j(27rfoto+0) 

e  (in  complex  notation),  minus  the  argument  of 

j2TT(f  o-f  jp,)to 

the  local  oscillator  phasor,  e  ,  minus 

2iTfjp,to  .  The  time  to  is  the  time  of  the  measurement, 
i.e.,  the  time  elapsed  since  the  chosen  origin  time,  t=0  , 
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when  phase  was  first  defined.  Clearly,  without  some  sort 
of  feedback  between  the  received  and  local  oscillator  fields 
which  defines  a  unique  t=0  ,  (as  is  done  in  a  PLL),  the 

measured  value  of  the  phase  will  be  a  function  of  what  t=0 
is  chosen.  However,  changes  in  the  phase,  whether  temporal 
or  spatial,  will  be  invariant  with  respect  to  the  choice  of 
the  origin  of  the  time  axis. 

Overview 

This  paper  begins  with  a  general  analysis  of  how  the 
heterodyne  or  intermediate  frequency  (IF)  signal  comes  about 
and  a  short  description  of  why  the  white,  gaussian  noise 
(WON)  assumption  is  valid  for  the  cases  of  interest.  Next, 
the  preprocessing  forced  on  the  system  by  the  CCD  detectors 
is  examined  and  a  Cramer-Rao  lower  bound  for  the  error  com¬ 
mitted  by  any  processor  utilizing  these  samples  to  estimate 
phase  is  derived.  Estimators  are  derived  for  various  as¬ 
sumptions  concerning  the  time  correlation  of  the  phase  pro¬ 
cess  using  both  maximum  likelihood  (ML)  and  state-variable 
approaches.  Finally,  the  various  processors  are  simulated 
and  the  results  are  presented  and  compared  to  the  perform¬ 
ance  of  a  first  order  linearized  phase-locked  loop  (PLL). 
This  is  a  standard  benchmark  for  comparisons  in  phase  esti¬ 
mation  problems.  Throughout  the  paper,  only  one  detector 
in  the  array  is  examined  and  no  spatial  information  is  in¬ 
corporated. 
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I I ,  Receiver  and  Preprocessor  Theoretical  Background 


In  this  chapter,  the  optical  heterodyne  receiver  is 
very  briefly  described  and  a  justification  for  the  WGN  as¬ 
sumption  is  presented.  Next,  the  preprocessor  structure  is 
examined  and  a  Cramer-Rao  (CR)  lower  bound  is  derived  for 
any  estimator  using  such  processed  samples.  Some  limiting 
assumptions  are  imposed  on  the  bandwidths  of  the  phase  pro¬ 
cess  and  the  amplitude  variations  of  the  IF  signal.  Finally, 
the  Cramer-Rao  bound  is  compared  to  that  of  a  processor 
utilizing  the  signal  without  integration  or  sampling. 

Optical  Receiver  Structure 

An  optical  heterodyne  reciever  is  shown  schematically 
in  Figure  1.  The  received  field  and  the  local  oscillator 
field  are  both  given  in  complex  notation  with  an  implied 
time  dependence  of  exp( -j2Tr  f  o  t  )  ,  where  fg  is  the  mean 

optical  frequency  of  the  received  field.  Each  element  of 
the  detector  array  responds  to  the  total  incident  optical 
power.  The  output  of  the  m,n-th  detector  is: 

"mn^^^  =  h-f^  /  f^R^^m’^n'^^  ^ 

"A  , 
d 

exp[-j2Trf  jpt]  I  Mx^dy^  (1) 
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Figure  1.  Receiver  Schematic  Diagram 

•  U£(x^.y^.t)exp[j2Trf  jj.t]  }]  dx^dy^  (2) 

^d 

*  I  •cos(2TTf^^t+9(x^,y„,t))}dx„dy„  (3) 


m’^n' 


Where  n  is  the  quantum  efficiency  of  the  detector,  h  is 
Planck's  constant,  U„  and  U,  are  the  received  and  local 
oscillator  (LO)  fields,  respectively,  measured  in  the  detec* 
tor  plane,  Up(x^,  y^  ,  t  )=  | Uj^(x^,  y^^ ,  t )  |  exp  [-j 0  (x^,  y^ ,  t )]  ,  and 
the  integration  is  performed  over  A^,  the  area  of  the  m,n-th 
detector.  The  local  oscillator  is  assumed  to  be  stabilized 
in  frequency  and  assumed  to  produce  a  uniform  plane  wave  so 
that  U.  (x„,y„t)=U  .  Assuming  that  0(x,y,t)  does  not  vary 

over  the  area  of  an  individual  detector,  factoring  out 


defining: 

"mn  “  HTtI  ‘IV=‘m->'n’‘>l“  +  l\l'’dx„<iy„ 

^  2lUj,(x^.y„.t)|.|Ujdx^dy„ 

^d 


reduces  equation  (3)  to: 

"mn'O  ■  "n,n[l"-''‘=°=<2’'lF*'"«mn<‘>>]  ‘®> 

Call  Y  the  modulation  index  and  i  the  amplitude.  It 

mn 

is  easily  seen  that  y^I  •  Since  only  one  detector  is  be¬ 
ing  examined  henceforth,  the  notation  will  be  simplified  by 

replacing  i  with  i„  and  a  (t)  with  e  (t)  .  Equa- 

mn  s  '^mn 

tion  (6)  is  the  signal  model  which  acts  as  the  input  to  the 
preprocessor  structure  ( integrate-and-dump  filters)  imposed 
by  the  CCD's. 

There  is  an  inferent  uncertainty  in  the  exact  arrival 
and  conversion  times  of  incident  photons  at  any  particular 
detector  imposed  by  the  statistical  nature  of  light.  These 
uncertainties  are  known  to  result  in  a  noise  process  associ¬ 
ated  with  the  signal  that  is  Poisson  distiibuted  (Ref. 5: 53). 
However,  if  the  arrival  times  of  the  photons  are  sufficiently 
close  together,  then  for  any  finite  bandwidth  detector,  the 
central  limit  theorem  forces  the  resulting  filtered  statis¬ 
tics  to  be  gaussian.  (Ref . 3 : 124 ) .  The  smaller  the  bandwidth 
for  a  given  arrival  rate,  the  more  closely  the  resulting  pro¬ 
cess  approximates  a  gaussian  distribution.  For  typical  local 
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oscillator  power  levels,  the  average  photon  arrival  rate  is 
so  high  that  even  relatively  wideband  detectors  will  exhibit 
gaussian  statistics  at  their  outputs  (Ref. 12).  For  the  low 
bandwidth  (compared  to  photoconductive  or  photovoltaic  de¬ 
vices)  CCD  devices  considered  in  this  paper,  the  gaussian 
assumption  is  extremely  good  for  local  oscillator  noise  lim¬ 
ited  cases. 

Preprocessor  Structure 

The  output  of  a  heterodyne  detector  with  noise  corrup¬ 
tion  is: 

r(t)  =  h(t)+n(t)  =  cos(2TTf  j^t  +  e  ( t ) )] +n(  t )  (7) 

where  h(t)  is  the  phase  modulated  carrier,  0(t)  is  the 
phase  process  to  be  estimated,  and  n(^)  is  a  white,  gaus¬ 
sian,  zero  mean  noise  process  which  accounts  for  the  detec¬ 
tor  noise.  Since  r(t)  is  the  output  of  a  total  power  de¬ 
tector,  it  clearly  can  never  be  negative.  But,  if  n(t)  is 
truely  gaussian,  then  there  is  a  finite  probability  that  it 
will  take  on  values  such  that  r(t)  does  indeed  become  neg¬ 
ative.  The  resolution  of  this  apparent  contradiction  lies  in 
the  fact  that  n(t)  is  never  precisely  gaussian.  The  noise 
n(t)  arises  as  a  result  of  the  fact  that  the  photon  arrival 
times  are  not  known  exactly,  even  if  the  field  is  determinis¬ 
tic  and  known.  For  a  field  with  intensity  (power  at  the  de¬ 
tector  plane)  I(t),  the  detector  output  will  be 
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r(t)  =  al(t)+n(t)  ,  where  a  is  a  constant  of  proportion¬ 
ality.  It  can  be  shown  that  the  signal  r(t)  is  a  Poisson 
process  with  rate  parameter  al(t)  (Ref. 5: 53).  The  mean 
of  the  process  is  given  by  oil(t)  and  so  is  the  variance. 
The  central  limit  theorem  requires  that,  with  a  finite 
bandwidth  detector,  as  the  number  of  photoelectron  conver¬ 
sions  per  second  approaches  infinity,  the  variable  (r(t)- 
al(t))//al(t)  will  become  gaussian  with  zero  mean  and  unit 
variance  (Ref. 3: 124).  That  is,  r(t)  becomes  gaussian  with 
mean  al(t)  and  variance  al(t)  ,  hence,  n(t)  becomes 
gaussian  with  mean  zero  and  variance  al(t).  However,  the 
photon  count  per  second  approaches  infinity  only  if  I(t) 
becomes  infinite.  This  would  force  the  mean  and  variance 
of  r(t)  to  infinity  as  well.  For  a  large  (but  finite) 
photon  arrival  rate,  the  process  will  becomes  very  close  to 
gaussian  around  the  mean,  but  cannot  be  accurately  modeled 
as  gaussian  far  into  the  "tails”.  This  is  so  since  a  Pois¬ 
son  process  is  defined  only  for  (discrete)  positive  argu¬ 
ments,  the  process  r(t)  ,  a  Poisson  process,  is  defined 
only  for  values  of  its  pdf  greater  than  zero.  That  is, 
f^^^ j (r(t ) )=0  for  r(t)<0  .  Therefore,  the  probability 

of  n(t)  taking  on  values  less  than  -al(t)  is  zero,  and 
n(t)  is  not  truely  gaussian.  Since  al(t)  is  very  large, 
the  truncation  of  the  tails  occurs  at  a  very  large  number 
of  standard  deviations  from  the  mean  for  even  moderate  in¬ 
tensities.  Hense,  the  assumption  that  n(t)  is  zero  mean 
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and  gaussian  is  very  good  for  large  fields  incident  on  the 
detector  (many  photo  electron  conversions  per  second  com¬ 
pared  to  the  bandwidth  of  the  detector).  See  Davenport  and 
Root  (Ref. 3: 124)  and  Karp  and  Gagliardi  (Ref. 5: 130). 

The  effect  of  using  a  CCD  detector  is  to  introduce  an 
integrate-and-dump  filtering  (IDF)  on  the  signal  before 
reaching  any  of  the  receiver  electronics  at  the  IF  stage. 

The  integration  period  of  the  CCD  is  user  controlable  as 
long  as  it  is  more  than  a  certain  minimum  value  set  by  the 
device  characteristics.  A  typical  minimum  value  might  be 
O.lysec  (Ref.2).  Designate  the  number  of  integration  per¬ 
iods  per  period  of  the  IF  by  .  The  output  of  the  IDF 

is  a  sequence  of  random  variables  with  the  form: 

r^  =  J  r(t)dt 

-  ^(i-l)+27rm(i-l) 

^(  i-^)+27rmi 

=  J  h(t)  +  n(t)  dt  (8) 

|-(i=|)+2TTm(i-l) 

where  T=l/fjp  .  The  factor  of  -T/2P^  occurs  in  the  in¬ 
tegration  limits  in  order  to  take  advantage  of  trigonometric 
simplifications  that  result.  The  factors  of  Strm  are 
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included  in  the  limits  to  indicate  that  it  is  not  necessary 
to  integrate  over  just  one  IF  period.  Since  the  IF  per¬ 
iod,  T  ,  is  typically  on  the  order  of  20  ns  or  less,  and 
the  CCD  connot  integrate  for  less  than  a  few  tenths  of  a 
nsec,  the  CCD  will  obviously  integrate  for  more  than  T/P^ 
seconds  per  sample.  Since  h(t)  is  2it-  periodic,  adding 
a  factor  of  2Trm  in  the  limits  does  not  change  the  samples 
except  that  the  variances  of  the  outputs  of  the  CCD's  are 
reduced  by  the  additional  integration  time  and  the  constant 

bias  (due  to  integrating  i  )  is  increased  deterministically. 

s 

Throughout  this  paper,  n  is  set  to  zero  simply  for  con¬ 
venience.  If  m^O  ,  the  restriction  on  the  upper  limit  to 
the  bandwidth  of  the  process  9(t)  ,  to  be  imposed  shortly, 

will  be  more  stringent.  As  shown  in  equation  (8),  the  se¬ 
quence  is  the  sum  of  a  signal  sequence  {h^}  and  a 

noise  sequence  {n^}  .  The  output  of  the  detector  is  ob¬ 

served  for  a  total  of  K  IF  periods,  the  length  of  the  ob¬ 
servation  interval  is  KT  second  and  the  length  of  the  ob¬ 
servation  vector  is  KP^  elements.  If  mj^O  ,  the  total  ob¬ 
servation  interval  is  K(mP^+l)T  .  The  observation  vector 
is  composed  of  K  subvectors  each  of  P^  elements: 


£  =  (ri ,r2 


^P^’^P^+1’  *  *  "^2P^'  •  •  *  '^(K-1)P^+1, 


•  • • » ^  ( 9 ) 
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(10) 


/  nT 

=  . Ik^ 

The  parameter  has  been  left  as  an  arbitrary  inte¬ 
ger.  If  there  is  a  value  of  that  minimizes  the  mean- 

squared  error  of  the  processor,  that  would  be  the  optimum 
value  of  P^  to  use  in  a  processor.  To  determine  if  such 
a  value  exists  and  what  that  value  is,  a  Cramer-Rao  (CR) 
lower  bound  is  computed. 

Cramer-Rao  Lower  Bound 

The  CR  lower  bound  for  the  unbiased  estimate  of  a  par¬ 
ameter  i,  is  given  by  (Ref.  17: 66); 

E[(5-0']^e'  =  -e[^  In  f(rU)]"*  (11) 

where  c  is  the  estimate  of  the  parameter,  r  is  the  total 
observation  vector,  f(r|^)  is  the  probability  density  func¬ 
tion  (pdf)  or  £  conditioned  on  knowing  the  true  value  c  , 
and  eQ  is  the  expectation  operator.  The  conditional  pdf 
f(r|0(t))  is  required.  The  signal  is  a  deterministic  func¬ 
tion  of  the  phase  and  the  noise  is  assumed  to  be  independent 
of  the  phase.  Since  the  noise  is  additive,  f(r|0(t))  is 
just  f(n)  shifted  to  a  mean  value  of  h 

Since  n(t)  is  white  and  gaussian,  and  the  integral 
is  a  linear  operator,  integrated  samples  of  n(t)  will  form 
a  gaussian  sequence.  It  is  easily  derived  that  {n^}  is  a 
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white  and  zero  mean  sequence,  as  well. 

E^n^]  =  E^  yn(t)dt]  =  jE^n(t)j  dt  =  0  (12) 


^  e[  Jn(t)dty  n(t')dt'j  (13) 

i  j 


=  J fE^n(t)xi(t')^dtdt'  =ff^&it-t')dtdt'  (14) 


NoT 


6.  . 


2Pt  "j 


(15) 


The  two-sided  spectral  density  of  n(t)  has  been  denoted 

No/2  .  The  appearance  of  the  Kroenecker  delta,  6. .  , 

^  J 

demonstrates  that  the  sequence  is  white,  as  expected.  The 
pdf  on  n  can  now  be  written  directly: 


f(n)  =  I  -  I  exp 

\/2iNot  / 


1 

2  ’NoT 


Similarly,  the  conditional  pdf  of  r|e(t)  is: 
2P.  \P^K 


f(r|0) 


\/2ttNoT/ 


exp 


2'NoT  S  (^i"*^i^‘ 
i=l 


(16) 


(17) 


From  equation  (11)  the  CR  bound  is: 

P. 


-Al/-  1  V*/-  V.  ^  2  xl  -  » 

=  -E  NoT  2^^*'i"^i^  ^ 

•  • 


(18) 
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I 

‘i 


( 


(19) 


Where  the  fact  that  ~  exploited  (since  n(t) 

is  zero  mean)  to  remove  the  term  involving  the  second  deri¬ 
vative  of  h^  .  Inserting  equation  (6)  into  equation  (8) 
gives : 


h.(e(t.  )) 


i  TY 
s 

27r 


cos(^(i-i)+0(t^))-cos(^(i-|)+0(t^)) 


(20) 


where  t^=T(i-J)/P^  .  It  was  also  assumed  in  this  integral 

that  9(t)  does  not  change  over  the  period  of  integration, 
else  the  integration  could  not  be  performed  without  explicit 
knowledge  of  the  functional  form  of  6(t)  ,  if  then.  Fur¬ 
thermore,  if  the  phase  changes  much  over  one  measurement  in¬ 
terval,  then  some  sort  of  average  of  the  phase  over  the  in¬ 
terval  will  be  estimated.  Consequently,  it  will  be  assumed 
that  the  bandwidth  of  the  process  9(t)  is  small  enough 
that  the  realization  being  examined  changes  very  little  over 
the  measurement  interval,  with  high  probability.  Equivalent¬ 
ly,  some  statistical  measure  of  the  bandwidth  of  the  phase 
process  is  assumed  known  a-priori  and  the  measurement  inter¬ 
val  is  chosen  so  that  the  phase  process  realization  will 
change  by  less  than  some  small,  predetermined  bound  during 
the  measurement  interval  with  some  high  degree  of  probabil¬ 
ity  chosen  by  the  designer.  If  the  m  in  equation  (8)  is 
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I 

I 

i 


! 


t 


1 


large,  this  can  be  a  rather  stringent  condition. 


With  these  assumptions,  the  CR  bound  becomes: 


=  P^rcos(x.)-cos(x.- 


X  =  ^  (i-1)  +  0 
t 


(22) 


Using  the  trigonometric  identity  for  cos(x)-cos(y )  and  apply¬ 
ing  the  results  of  Appendix  A  to  the  summation  results  in; 


=  i^(|^)^[2P,^-K.sin^(^)]-' 

S  L  ^  J 

_  pNo  2  1 

[2TK  *  (igY)"J  *  sin^(TT/P^) 


N  2 

As  P^'*’”  ,  this  equation  reduces  to  ^^~2lgT‘  ( i  y')'^'  which 

s 

is  the  same  as  the  error  for  the  linearized  first  order  PLL 
(Ref. 18: 93),  or  the  CR  bound  for  a  system  which  has  access 
to  the  non-integrated,  non-sampled  waveform.  For  P^  fin¬ 
ite,  the  error  bound  will  be  greater  than  that  for  the  PLL 
by  a  factor  of  ^^/P^  sin^(^)  .  In  all  simulations,  P^= 

4  simply  to  reduce  computer  time  requirements.  For  larger 
P  ,  an  appropriate  scale  factor  can  be  introduced  into 
the  simulated  errors  by  using  equation  (24). 
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III.  Processor  Determination : 
Maximum  Likelihood  Estimation 


Three  specific  cases  of  the  estimation  problem  can  be 
identified; 

(i)  0(t)  is  a  constant,  0o  ,  for  all  time.  It  may 
be  construed  as  a  single  realization  of  a  random  variable 
uniformly  distributed  between  -n  and  tt  radians. 

(ii)  0(t)  is  a  random  process  with  unknown  statistics 

(iii)  0(t)  is  a  random  process  with  known  statistics. 
The  first  two  cases  are  analyzed  in  this  chapter  and  it  is 
shown  that  the  optimum  processor  is  a  discretized  phase-lock 
ed  loop.  Because  of  the  nature  of  the  integrate-and-dump 
preprocessing  imposed  on  the  measurements,  it  may  not  be  ob¬ 
vious  that  the  optimum  ML  processor  still  has  a  PLL  struc¬ 
ture.  The  third  case  is  analyzed  in  the  next  chapter. 


Measurement  Model  Equations 

Before  continuing,  it  is  convenient  to  derive  the  exact 
form  of  h^(0)  and  its  first  two  derivatives  with  respect 
to  0  .  Using  equation  (6)  in  equation  (8)  gives  h^(0)  , 

which  can  be  differentiated  directly. 


h^(d) 


igT  ^  ^  sin(J-)cos(|^(i-l)+0) 

m  • 


15 


6h  (6)  i  TY  ^  2tt 

^  sin(pl)sin(|-(i-l).e) 


6^h  (6)  i  Ty  ^  21, 

-TF-  =  -  -f-  sin(^)cos(|l(i-l)-He) 


These  equations  are  used  many  times  in  the  derivations  that 
follow  and  will  be  referred  to  quite  often.  These  equations 
are  needed  because  the  first  derivative  is  always  necessary 
in  the  derivation  of  the  ML  processor.  The  second  derivative 
will  be  used  when  the  Kalman  filter  is  derived  because  the 
nonlinear  measurements  will  be  linearized  in  a  Taylor  series 
expansion  and  the  second  derivative  appears  in  the  deriva¬ 
tion  . 


Maximum  Likelihood  Approach:  Case  (i) 

The  ML  estimate  of  9o  is  the  value  of  9n  that  maxi¬ 


mizes  the  conditional  density  f(rj6o)  .  Since  f(£l9o) 


is  gaussian,  the  ML  estimate  is  the  value  of  9 


0  >  o  0  > 


that  solves: 


In  f(£|9o)  -  0 


Taking  the  logarithm  of  equation  (17)  and  differentiating 


gives : 


V'  °"i 
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Inserting  equations  (26)  and  (27)  results,  after  some  trig¬ 
onometric  identities  are  applied: 


0  =  sin(|-(i-l)+0o )  -  ^ 


i  T 

sin(^)sin(|^-(i-l)  +  f 
t  t  t 


Yi  T 
s 

2Tr 


sin(^) 

t 


sin(ii(i-l)+20o) 

t 


(30) 


The  trigonometric  terms  in  the  last  two  summations  are 
shown  to  sum  to  zero  in  Appendix  A.  The  estimator  reduces 
to : 

KPt 

Y.  r.  sin(|l(i-l)+0)  =  0  (31) 

i=l  ^ 

This  processor  is  a  discrete  analog  to  the  familiar 
PLL  and  can  be  implemented  in  closed  loop  form  for  any  val¬ 
ues  of  K  and  P^  .  A  block  diagram  is  shown  in  Figure 
2,  The  dotted  box  encloses  a  discrete  IDF(DIDF). 

Maximum  Likelihood  Approach:  Case  (ii) 

Suppose  that  0(t)  is  not  constant  over  the  entire  ob¬ 
servation  interval.  In  light  of  Appendi.x  A,  let  the  phase 
be  approximately  constant  over  one  period,  T  ,  but  vary  from 
period  to  period  in  an  unknown  manner.  Equation  (31)  be¬ 
comes  : 
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If  equation  (32)  is  to  be  true  for  any  value  of  K  ,  then 
the  inner  sum  must  be  zero  for  every  m  .  The  estimator 


is  then; 


E  '■i  +  -  0 

i=l  ^ 


The  estimator  produces  a  sequence  of  estimates, 

might  be  expected,  0^  is  estimated  independently  of  any 

/N 

past  value  of  0  or  0  and  only  data  from  the  current 
mm 

measurement  interval  ^mT,  (m+l)Tj  is  used  to  perform  the 

T 

estimation.  That  is,  r=(r  ,r  , . , , , rp^ )  .  This  is  direct¬ 

ly  a  result  of  the  fact  that  {0i,}  is  (or  is  assumed  to  be) 


an  uncorrelated  sequence; 


0  kl  • 
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IV.  Processor  Determinat ion :  State  Variable 


Analysis  and  Kalman  Filtering 


If  some  information  concerning  the  time  correlation  of 
the  process  6(t)  is  already  known,  it  should  be  possible 
to  obtain  a  better  estimate  (lower  variance  estimate)  of 
the  process  by  incorporating  the  a-priori  information  into 
the  processor.  In  this  chapter,  state  variable  models  and 
nonlinear  Kalman  filtering  are  applied  to  the  problem  of 
incorporating  a-priori  information  into  the  processing  al¬ 
gorithm.  A  few  words  are  said  about  state  variable  equa¬ 
tions  in  general  and  then  the  arguments  are  specialized 
for  the  measurement  model  of  interest. 

State  Variable  Equations 

The  phase  process  is  assumed  to  be  the  output  of  a  lin¬ 
ear,  time  invariant  system  driven  by  white  gaussian  noise. 
This  system  is  described  by  n  state  variables  in  a  vector 
differential  equation  (Ref . 13 : 556 ) . 


x(t)  =  F  x(t)  +  Gu(t) 

(34) 

9(t)  =  c”^  x(t) 

(35) 

c”^  =  (0,0,. ..,0,1) 

(36) 

where : 
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x(  L. )  n-dimensional  state  vector 
u(t)  m-dimensional  noise  vector 
F  nxn-dimensional  system  matrix 

G  nxm-dimensional  noise  system  matrix 

The  solution  to  equation  (34)  is  (Ref . 13 : 556) : 


x(t)  =  <))(t-t  0  )x(t  0  )  +  /  ( t-t  ^  )Gu(t  '  )dt  ■ 


Where  (t)(t-to)  is  called  the  system  transition  matrix  and 
solves  equation  (38)  with  initial  condition  ( 39) (Ref . 13 : 557 ) 


(Kt-to  )  =  F(t)(t-to ) 


<t)(0)  =  I 


Since  the  data  to  be  processed  is  discrete,  a  discrete  for¬ 
mulation  of  equation  (37)  is: 


x(kT)  =  Xk  =  [<t>(T)]^  xo  +  S  [<J>(T)]^  ^G  u^  (40) 


This  equation  can  be  put  into  a  recursive  form  with  only  mod¬ 
erate  algebraic  manipulation  (Ref . 13 : 557) : 


^k-1 


For  the  time  invariant  case,  the  solution  to  equation  (38)  is 
exp(F( t-t  0 )  )=(|)(  t-t  0 )  which  becomes,  for  discrete  measure¬ 
ments  (Ref . 13:558) : 
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<p  =  (t)(T)  =  exp^Fxj 


(42) 


It  will  now  be  assumed  that  the  phase  does  not  change 
during  an  IF  period  (or  changes  very  little)  so  that  the 
phase  process  is  described  by  a  sequence  of  random  variables 
{6j^}  or,  alternatively,  >  where  the  variables  are  sep¬ 

arated  by  T  seconds.  This  will  allow  taking  advantage  of 
the  trigonometric  simplifications  proven  in  Appendix  A. 

-k  ^  -k  ""  £k 

Where  nj^  is  a  white,  gaussian  vector  with  statistics: 


The  dimensionality  of  r,  h,  and  n  is  .  The  matrix 

I  is  a  identity  matrix.  The  sequence 

white  and  gaussian  with  covariance  matrix  Q  .  The  sequences 
{nj^}  and  a-re  independent.  These  assumptions,  plus  the 

initial  conditions  are: 


E 


®[Hk]  ■  2 

E[uk3ll]  -  <i‘kl 

(45) 

r  Ti  „ 

Ti 

E[xkHl]  = 

(46) 

e[(x-xo  )(x-Xo  )^  =  Po 

(47) 
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Given  the  received  sequence  .  an  "optimum"  es- 

timate  of  the  state  is  desired  from  which  an  optimum 

estimate  of  the  phase  can  be  derived,  or  any 

other  element  of  the  state  vector.  Unfortunately,  no  gen¬ 
eral  solution  to  this  problem  exists  for  arbitrary  nonlin¬ 
ear  functions  h( ♦ )  in  either  the  continuous  or  discrete 
case.  Exact  solutions  do  exist  for  linear  h( • )  and  ap¬ 
proximate  solutions  do  exist  for  a  nonlinear  ( •  )  which 
has  been  linearized  by  a  Taylor  series  expansion  about  the 

/V 

best  estimate  of  the  state,  Xj^  .  Sage  solves  the  prob¬ 
lem  using  a  discrete  invariant  embedding  procedure  (Ref. 14: 
450)  to  give  the  equations  for  a  maximum  a-posteriori  (MAP) 
estimate : 

Ik  " 


ik+l  ■  '**k  •’k*!' 


-T  6 


6Xk 


h  (4k)\’[ik+i-ii(<^>ik>] 


k+1 


’  -T  ^-T 

6 

<P  -<t> 

6''T 

^k_ 

-^h’^(<t,£  )R  ' 


r^^l-h(4,Xj^) 


^■'Pk 

(50)- 


^k  "  GQg'^4."'^  (51) 

When  h(4>Xj^)"ll(£  »  this  can  reduce  to  (see  Appen- 

dix  B) : 

-k  "  Ilk 
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where : 


CNR 


2  /  2T 


(62) 


i  Y 
s 


sin(^)cos(|i(i-l)  +  0j^^l)+nj^^l(i) 


(63) 


2  =  No 
n'  2P^T 


(64) 


The  carrier-to-noise  ratio  (CNR)  thus  defined  is  the  ratio 
of  the  signal  power  to  the  noise  power  in  the  total  band¬ 
width  used  in  producing  the  estimate  (1/T)  .  This  is  the 

same  CNR  used  by  Viterbi  (Ref. 18: 93)  and  thus  allows  direct 
comparisons  with  his  results  for  the  PLL,  In  addition,  the 
CNR  gives  an  intuitive  feel  for  how  "good"  the  measurements 
are,  or  how  badly  they  are  corrupted  by  measurement  noise. 


Processor  Ambiguity 

Appendix  B  shows  that  the  processor  equations  can  be 
further  simplified  to: 


-k+1 


\+l  "ks] 


(65) 


(66) 
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where : 


6  =  :^-sin^  (^)-CNR 


(67) 


var 


=  var 


=  6 


(68) 


The  sequences  also  zero  mean.  For 

any  plant  model  which  is  one  dimensional,  — k~^k”^k 
T 

cc  =1  .  For  this  case,  plugging  equation  (66)  into  (65), 

assuming  the  noise  is  negligible,  and  dividing  top  and  bot- 
tom  by  cos(9j^^j^-9^)  : 


®k+l  =  ®k  - 


tanOk+l-ep 


6  ■  ^  (})Pj^(t) '^+Q  "  )  "sec  ( 0  e  ' )  + 1 


(69) 


For  large  CNR,  the  estimate  should  be  quite  good  and,  with- 

/s 

in  modulo-2Tr,  ®k+l~®k~*^  *  sec(0)-l  ,  equation  (69) 

reduces  to,  approximately: 


®k+l  =  ®k  - 


-I 


1+8  ((})Pj^4)^+Q") 


- 1 


(70) 


However,  the  function  tan(»)  exhibits  a  modulo-ir  ambiguity, 
so  that  tan(0j^^j^-(§'  +  PTr))=tan(ej^^j^-e')  ,  where  p  is  any 
Integer.  Thus,  the  estimator  will  be  unable  to  determine 
the  phase  better  than  modulo-ir 

The  above  arguments  were  for  a  one-dimensional  state 
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vector  with  high  CNR.  Simulations  have  confirmed  this  re¬ 
sult  for  smaller  CNR's  as  well. 


Constant  Phase 

The  case  of  a  constant  phase  plant  is  the  simplest 
case  to  model.  The  state  vector  has  but  one  component  and 
the  matrices  collapse  to  scalars.  The  state  variable  equa¬ 
tions  become; 


9 


k+1 


(71) 


(t)  =  l  G  =  0  R  =  Q  =  arbitrary  (72) 

^^t 

The  estimator  or  Kalman  filter  equations  are,  from  (65)  and 

(66); 


'k+1 


=  9,.  -  P 


k+1 


[6-sin(e-Gj^)+nJ 


(73) 


Figure  (3b)  is  a  block  diagram  of  this  processor.  Figure 
(3a)  shows  the  processor  in  the  form  given  by  equations  (57) 
and  (58). 

Random  Phase  Walk 

The  random  walk  is  modelled  by  the  state  variable  equa¬ 
tion  ; 

®k+l  =  °k  ^  G^k 

27 
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Figure  3a.  Constant  Phase  Processor  Structure  (Realization  1). 
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]=Q=l  , 

of  the  noise  term  is 
$.=1  ,  and  R=NoT/2P^ 

filter  equations  are: 


Let 


var 


h 


for  simplicity,  so  that  the  variance 
reflected  in  the  value  of  G  .  Also, 
as  in  the  constant  phase  case.  The 


k+1 


=  0,  -  P 


k+ 


B-sin(0k^l-0j;)+nj^J 


(76) 


(VG')  '  ^^°^(®k+l-®k>^"kc 


(77) 


Figure  (4b)  shows  this  processor  while  Figure  (4a)  shows 
the  implementation  according  to  equations  (57)  and  (58). 
Both  forms  of  the  filter  are  shown  because  while  the  form 
of  equations  (65)  and  (66)  gives  a  good  intuitive  feel  for 
how  the  filters  process  the  estimation  error,  the  form  of 
equations  (57)  and  (58)  is  the  actual  form  a  hardware  pro¬ 
cessor  would  take.  This  is  true  because  actual  phase  is 
quite  clearly  unavailable  to  the  processor. 


Demodulation  of  FM  Signals 

Suppose  that  the  information  of  interest  in  a  state 

T 

vector  is  given  by  m(t)=c  x(t)  .  Suppose  also  that  the 

phase  of  the  carrier  is  given  not  by  m(t)  ,  but  by  0(t)= 

^  .  T 

X.^m(T)dT  +  0o  so  that  0(t)=m(t)=Ac  x(t)  where  X  is  called 

the  frequency  modulation  index.  A  new  state  vector  can  be 

constructed  from  the  old  state  vector  and  0(t)  .  The  new 

state  variable  equation  is: 
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X  (t) 


x(t) 

'f  0  “ 

'x(t) 

"g  0  ‘ 

u(t ) 

= 

+ 

_o(t)_ 

1 - 

lo 

H 

o 

1 _ 

0(t)_ 

Ug(t)_ 

x^(t)  =  F'  x'(t)  +  u'(t)  (79) 

Where  F_  has  been  added  to  the  F'  matrix  to  account  for 

U 

possible  variations  in  6(t)  not  due  to  the  FMing  of  the 
carrier  by  m(t).  In  the  discrete  formulation: 

^k+1  =  +  'J)"G"u^  (80) 


(81) 


The  plant  equation,  and  hence  the  filter  equation,  is  at 
least  two  dimensional,  but  this  is  no  theoretical  compli¬ 
cation. 


A  practical  case  is  when  the  frequency  is  a  constant 
over  the  observation  interval.  This  could  come  up  when 
estimating  phase  in  the  presence  of  an  unknown  doppler 
shift.  Model  the  system  with  two  state  variables  and  sup¬ 
pose  the  phase,  in  the  absence  of  the  doppler  shift,  would 
be  executing  a  random  walk.  Then: 


1 

o 

o 

x(t ) 

1 

o 

o 

_ 1 

'u(t  ) 

+ 

1 

o 

_ 1 

1 

o 

_ 1 

1 — 

o 

o 

o 

1 _ 

.UQ(t). 

(82) 
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where ; 

It  is  now  possible  to  estimate  both  phase  and  frequency. 

The  doppler  shift  is  given  by  ^k+l~—  -k+1  phase  is 

given  by  ’  where  d^=(l,0) 

Phase  Estimation  with  Frequency  Uncertainty 

Imagine  a  case  where  it  is  desired  to  estimate  the  phase 
of  an  optical  field,  but  the  local  oscillator  is  known  to  be 
unstable.  Alternatively,  the  LO  may  be  stable  and  the  re¬ 
ceived  field  frequency  fixed,  but  the  integration  limits  may 

not  be  controlled  precisely.  The  uncertainty  in  the  limit 
T 

times,  TT-Ci-a)  ,  can  be  modelled  as  an  uncertainty  in  fT^^^ 

Ir 

1 
T 
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The  frequency  uncertainty  can  be  incorporated  into  the  model 
of  equation  (80)  by  simply  allowing  GfO  .  Unfortunately, 
the  state-variable  models  assume  a  driving  force  which  is 
white,  and  the  frequency  uncertainty  of  a  local  oscillator 
is  not  even  nearly  white  for  any  case  of  interest.  The 
state-variable  models  can  accomodate  non-white  driving  func¬ 
tions  by  driving  a  "coloring  filter"  with  white  noise  and 
allowing  the  output  of  the  coloring  filter  to  drive  the  sys¬ 
tem  of  interest.  The  state  variable  equations  for  the  col¬ 
oring  filter  are  then  used  to  augment  the  original  state 
variable  equation  (Ref. 7:180).  The  coloring  filter  is  de¬ 
scribed  by; 

x^(t)  =  F^x^(t)  +  G^u^(t)  (88) 

u(t)  =  HjX^(t)  (89) 

Where  the  notation  in  equation  (89)  is  meant  to  indicate 
that  only  the  u(t)  portion  of  u^(t)  in  equation  (79) 
is  augmented  while  UQ(t)  (and,  for  that  matter,  any  por- 
tions  of  u(t)  that  are  already  white)  is  not  an  output 
of  the  system  augmented  system  is  described 

by; 
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(90) 


(91) 


The 


x(t)=  F  •x(t)+G 
—a  a  —a 

discretized  equation  is 


•  u  ( t ) 
a  —a 


given  by: 


X  (k+1)  = 

— cl 


i  X  (k)+4)  G  u  (k) 


(92  ) 


$a  =  exp[Fj]  (92) 

Where,  for  ease  of  reading,  the  subscript  k  has  been  re¬ 
placed  by  the  argument  k  .  The  filter  ecjualioiis  are  sim¬ 
ply  equations  (65)  and  (66),  or  equations  (57)  and  (58), 
with  all  quantities  replaced  by  their  associated  augmented 
counterparts.  Note  carefully  that  u  (t)  is  now  a  white 
process  driving  the  augmented  system.  This  is  as  far  as  it 
is  possible  to  go  without  assuming  certain  correlations 
for  Xf(^) 

It  is  cautioned  that  it  was  assumed  that  e(t)  could 
be  modelled  as  the  output  of  a  system  driven  by  white  noise. 
If  this  is  not  the  case,  then  equation  (89)  is  modified  so 
that  u  (t)=H^Xj(t)  .  Also,  it  was  assumed  that  u(t)  is 

the  output  of  a  system  all  of  whose  inputs  were  non-white. 

If  this  is  not  the  case,  then  some  elements  of  u(t)  are 
not  included  in  (89),  are  not  augmented,  and  enter  equation 
(90)  in  their  original  form.  A  further  caution  is  that 
since  the  white  noise  driving  function,  Uj(t)  ,  is  gaus- 
sian,  the  resulting  density  for  u(t)  ,  the  frequency 


36 


uncertainty,  is  also  gaussian  for  any  linear  plant  model. 
This  is  because  a  gaussian  process  remains  gaussian  after 
being  transformed  by  a  linear  system.  So,  if  a  white,  gaus¬ 
sian  process  u^Ct)  drives  a  linear,  time  invariant  system 
to  produce  an  output  u(t)  and  u(t)  drives  a  laser  mod¬ 
eled  as  a  voltage  controlled  oscillator  centered  at  fo  , 
the  pdf  of  the  laser  output  frequency  is  the  same  as  the  pdf 
of  u(t)  ,  i.e.,  gaussian.  A  nonlinear  plant  would  be  re¬ 
quired  in  order  to  accurately  model  a  homogeneously  broad¬ 
ened  laser  line,  which  is  Lorentzian,  since  only  through  a 
nonlinear  transformation  can  a  gaussian  process  UjCt)  be 
transformed  into  a  nongaussian  process  u(t)  .  A  Lorentzian 
random  variable  is  also  known  as  a  Cauchy  random  variable. 
The  pdf  of  a  Cauchy  distributed  random  variable  is  given  by 
(Ref . 19:48) ; 


f 


X 


_ a 

a'^  +  Cx-m  )'^ 
'  x' 


(94) 


(Technically,  of  course,  none  of  the  moments  of  this  density 
exist.  That  is,  /x’^f^(x)dx  does  not  converge  over  infin¬ 
ite  limits  for  n>0  ) .  Of  course,  an  approximation  would 
still  be  made  later,  when  the  estimator  equations  are  lin¬ 
earized  in  a  Taylor  series  expansion  and  this  will  result 
in  a  different  modeling  error. 
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V.  Simulation  of  Kalman  Ki Iter  Processor 

In  this  chapter  the  signal  models  of  chapter  IV  are 
simulated  and  the  results  of  the  simulations  are  discussed. 
Comparisons  are  made  with  CR  bounds  and  a  linearized  first 
order  PLL.  Some  interesting  behavior  on  the  part  of  the 
processor  is  noted  and  explained. 

Simulation  Parameters 

A  subroutine  was  written  in  Fortran  IV  which  performs 
a  Monte  Carlo  simulation  of  the  processor  of  interest.  The 
simulator  performs  an  ensemble  average  over  a  user  deter¬ 
mined  number  of  realizations  of  the  measurement  noise  pro¬ 
cess.  The  subroutine  is  listed  in  Appendix  C. 

For  the  case  where  the  IF  is  fixed  and  known  (no  dop- 
pler  shift  is  present),  the  phase  was  modelled  in  one  di¬ 
mension  and  two  models  were  investigated.  The  first  model 
was  for  the  phase  a  constant.  The  second  model  was  for 
the  phase  executing  a  random  walk  around  the  unit  circle. 
These  two  models  were  chosen  because  they  represent  two  ex¬ 
tremes  in  the  degree  of  correlation  of  the  phase  process 
and  are  typically  the  models  encountered  in  the  literature 
(Ref. 13:939,  Ref. 8:48).  Neither  the  constant  phase  plant 
model  nor  the  random  walk  phase  plant  model  may  be  adequate 
for  any  particular  case,  but  since  they  bracket  the  process 
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correlations  likely  to  be  encountered,  they  can  serve  as 
upper  and  lower  bounds,  respectively,  on  how  well  the  fil¬ 
ter  will  perform  in  practice. 

Also  simulated  was  the  case  where  the  phase  is  execut¬ 
ing  a  random  walk  and  there  is  a  constant,  unknown,  doppler 
shift.  The  errors  in  jointly  estimating  both  phase  and  fre¬ 
quency  are  investigated. 

In  all  simulations  100  sequential  time  estimates  of 
the  phase  process  are  computed  and  each  estimate  is  en¬ 
semble  averaged  over  100  realizations  of  the  noise  process. 
One  hundred  realizations  of  the  noise  is  not  very  many,  but 
it  is  common  practice  in  Monte  Carlo  simulations  to  use  a 
number  of  realizations  that  is  relatively  small  compared  to 
what  the  central  limit  theorem  might  deem  prudent  for  a 
reasonable  confidence  in  the  statistics  (Ref. 16:939,  15:382, 
10.67).  Compared  to  the  references  cited,  100  is  a  rather 
large  number  of  realizations.  The  reason  that  so  few  real¬ 
izations  are  used  is  the  quantity  of  computer  time  consumed. 
One  run  of  the  subroutine  with  100  realizations  and  100  re¬ 
cursions  will  consume  about  10  to  15  seconds  of  CP  time. 

For  a  0.9  probability  that  the  phase  error  estimate  is  with¬ 
in  0.1  radians  of  the  actual  phase  error,  the  central  limit 
theorem  would  require  809  realizations.  The  weak  law  of 
large  numbers  would  require  3290  realizations  (Ref. 10:67). 
The  Chernoff  bound  requires  1515  realizations  (Ref. 19:97). 
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Constant  Phase  Processor 

For  the  constant  phase  processor,  the  phase  was  mod¬ 
elled  as  a  random  variable  uniformly  distributed  between 
-  Tf  and  77  radians.  Initial  values  were  then  Xo=^o=0  , 
Po=tt^/3.  Po  is  the  variance  of  the  a-priori  estimate  0o  . 
The  state  vector  is  one-dimensional.  A  value  of 
integrals  per  IF  period  T  was  used. 

Figure  5  shows  the  simulation  results  for  an  actual 
angle  of  6  =0  and  three  different  CNEE'S  of  +30db,  +10db, 
and  -lOdb.  For  CNR=30db,  the  processor  converges  almost 
immediately  to  the  correct  value.  For  CNR=-10db  ,  the 
processor  fails  to  converge  at  all  after  100  recursions. 
Because  the  derivation  of  the  Kalman  filter  performed  by 
Sage  discarded  all  terms  in  a  Taylor  series  expansion  of 

/s 


best  estimate  of  before  measurement 


Xj^  and  Pj^  beyond  the  first  order,  the  estimator  does  not 
perform  well  in  a  high  noise  environment.  In  a  high  noise 
environment,  the  assumption  of  approximate  linearity  of  the 
processor  around  a  best  a-priori  estimate  trajectory  (the 

is  incorp¬ 
orated,  ^  violated  and  consequently  the  filter 

no  longer  adequately  extracts  information  from  the  measure¬ 
ments.  As  a  result,  the  processor  fails  to  converge.  In 
the  low  CNR  environment,  the  processor  variance  does  not 
settle  to  a  steady  state  value  as  it  does  in  the  high  CNR 
environment,  but  rather  it  oscillates  about,  never  reaching 
a  steady  state.  The  processor  simply  tracks  the  noise. 
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For  CNH=10db,  the  processor  docs  converge  after  about  25 
recursions,  but  it  converges  to  the  wrong  value,  about 
5ioy=-0.01  radians  (-0.57®).  The  processor  is  biased. 

Many  texts  either  assume  that  nonlinear  Kalman  filters  are 
unbiased  or  state  without  proof,  justification  or  condi¬ 
tioning  that  the  bias  is  zero  (Ref . 14 :433) .  As  discussed 
in  Maybeck  (Ref. 8),  these  claims  are  false.  Maybeck  dis¬ 
cusses  methods  for  reducing  the  bias  of  the  processor  by 
using  second  order  approximation  filters  (retaining  Taylor 
series  terms  up  to  order  2),  or  by  using  first  order  filters 
(like  this  one)  but  adding  bias  correction  terms  taken  from 
the  second  order  filter  equations.  These  methods  result  in 
somewhat  more  complicated  filter  algorithms  and  were  not 
used  in  this  research,  however,  processor  biases  will  be 
pointed  out  as  they  are  encountered. 

To  investigate  the  bias  problem  further,  Figure  6  plots 
the  estimation  error  after  100  recursions,  ®ioO~®  ’ 
should  be  very  close  to  the  actual  bias,  for  three  CNR's. 
First,  notice  that  for  CNR=-10  db,  the  bias  varies  linearly 
(approximately)  with  the  actual  angle  0  .  This  is  because 

the  processor  is  tracking  the  noise  which  is  zero  mean  and 
tends  to  guess  0=C  ,  the  a-priori  estimate,  on  the  average 

Also  notice  that  even  for  a  very  high  CNR,  30db,  the  process 
or  guesses  0=0  ,  on  the  average,  whenever  0=±n/2  .  This 
phenomenon  is  a  direct  result  of  the  modulo-ir  ambiguity  of 
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the  processor.  If  the  actual  angle  is  +Tr/2  ,  the  process¬ 

or  may  determine  an  estimate  that,  due  to  noise  corruption, 
is  0± e  ,  where  c>0  is  the  error.  But  ^+e(or  -^-e)  is 
outside  the  ambiguity  range  of  the  processor  and  so  it 

X  TT 

guesses  il''®  processor  "wraps 

around"  the  semicircle  and  the  average  over  many  realizations 

is  zero.  This  explains  why,  even  though  the  estimate  seems 

poor  from  the  examination  of  a  sample  mean,  the  variance,  P,  , 

iC 

still  gets  small  as  k  gets  large.  It  also  explains  why  the 

bias  rises,  for  CNR=10db,  as  0-*^  .  With  the  higher  noise 

level,  as  the  angle  gets  closer  to  -  ^  ,  the  processor  will 

"wrap  around"  more  often  and  the  bias  will  rise.  Similarly, 

for  a  fixed  angle,  the  processor  will  "wrap  around"  more  often 

for  a  smaller  CNR,  so  the  bias  increases  as  CNR  decreases. 

_  1  _2 

Figure  7  is  a  plot  of  the  inverse  variance  ^ 

versus  the  CNR,  Also  plotted  is  the  theoretical  CR  bound^^*^ 
and  the  variance  for  the  linearized  PLL.  The  processor  sig¬ 
nificantly  underperforms  the  PLL,  although  this  gap  would 
narrow  if  P^  was  increased.  The  processor  almost  reaches 
its  CR  bound  with  equality,  however.  After  100  recursions, 
the  estimator  is  "almost  efficient". 

Random  Phase  Walk 

The  next  model  simulated  is  for  the  phase  executing  a 
random  v/alk  around  the  unit  circle.  For  this  case,  the 
variance  of  the  walk,  or  diffusion,  parameter  was  set  at 
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igure  7.  Constant  Phase  Processor  Variances 


Q=0.005  so  that  tlie  size  of  the  "steps"  in  the  walk  would 
be  small  enough  that  the  phase  can  be  assumed  approximate¬ 
ly  constant  over  each  IF  period.  For  this  value  of  Q  , 
the  magnitude  of  less  than  0.1  radians  (5.7°) 

with  probability  0.84.  Only  one  realization  of  the  phase 
process  was  examined  and  the  processor's  performance  was, 
again,  averaged  over  100  noise  process  realizations.  Ideal¬ 
ly,  the  processor's  performance  should  be  averaged  over 
many  realizations  of  the  phase  process  as  well  as  averag¬ 
ing  over  many  realizations  of  the  noise  process.  This 
averaging  was  not  performed  because  of  the  amount  of  compu¬ 
ter  time  that  it  would  have  consumed. 

Figure  8  shows  the  error  performance  of  the  processor 
plotted  for  100  recursions  at  3  CNIi's.  The  quantity  plot- 
ted  is  the  estimator  error,  •  The  reduction  in  per¬ 

formance  at  progressively  lower  CNR's  is  self  evident.  Note, 
however,  that  the  processor,  at  lOdb  CNR,  tends  to  made  more 
negative  errors  than  positive  errors.  This  is  another  in¬ 
dication  of  the  processor  bias.  The  bias  is  even  more  evi¬ 
dent  for  the  CNR=-10db  case.  For  this  particular  realiza¬ 
tion  of  the  walk,  the  phase  managed  to  walk  as  far  as  +1.4 
radians,  which  is  quite  close  to  .  For  a  walk  realiza¬ 
tion  which  does  not  walk  as  far  from  0=0  ,  the  error,  on 

the  average,  would  be  expected  to  be  less.  Simulations  con¬ 
firm  that  this  is  true,  although  the  processor  still  makes 
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significantly  more  errors  for  the  lower  CNR. 

Figure  9  plots  the  inverse  error  variance  versus  the 
CNR  along  with  the  CR  bound  and  the  theoretical  results 
for  the  linearized  PLL.  The  scale  on  the  vertical  axis 
has  been  changed  relative  to  Figure  7.  This  change  re¬ 
flects  the  increased  bandwidth  of  the  process  and  the  cor¬ 
responding  increase  in  the  bandwidth  of  the  PLL  loop  fil¬ 
ter.  The  loop  filter  time  constant  has  gone  from  KT,  where 
K=k  =100,  to  T  .  The  fact  that  the  processor  seems  to 

outperform  its  CR  bound  can  be  explained  on  the  basis  of 
the  fact  that  the  computation  of  the  CR  bound  in  equation 
(24)  set  K=1  .  That  is,  the  CR  bound  is  for  a  processor 
which  does  not  take  into  account  the  a-priori  information 
on  the  corre]ation  of  the  phase  process.  The  processor 
simulated  here  does  use  this  a-priori  information  and  it 
gives  a  lower  variance  estimate  as  a  result. 

Since  the  processor  is  biased,  and  the  CR  bound  was 
derived  for  an  unbiased  estimator,  it  might  be  argued  that 
the  processor  outperforms  its  CR  bound  because  the  bound 
was  derived  for  the  unbiased  case.  The  CR  bound  for  a  bias 
ed  estimator  is  given  by  multiplying  the  right  hand  side  of 
equation  (11),  (the  bound  for  an  unbiased  estimator)  by 
( l-db( 0 ) /do ) ^  ,  where  b(0)  is  the  processor  bias  as  a 

function  of  0  ( Ref . 17 : 146 ) .  From  Figure  6,  for  moderate 

to  high  CNR's,  db(0)/d6=O  for  all  angles  not  close  to  ± 
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At  low  CNR's,  the  Kalman  filter  also  outperforms  the  PLL. 
This  is  also  due  to  the  fact  that  the  PLL  does  not  incorp¬ 
orate  a-priori  information  about  the  plant  dynamics,  where¬ 
as  the  Kalman  filter  does  incorporate  this  information. 


Phase  Estimation  with  Doppler  Shift 

It  was  mentioned  earlier  that  phase  and  frequency  can 
be  estimated  jointly.  The  state  vector  becomes  two  dimen¬ 
sional  with  components  Xj=doppler  shift  and  X2=phase.  The 
phase  process  was  modelled  as  a  random  walk  with  diffusion 
parameter  Q=0.005  .  A  static  model  of  the  doppler  shift 

was  assumed  with  a  shift  of  +2MHz.  The  IF  was  chosen  at 
l/T=50MHz.  If  the  filter  is  in  a  1.06y  Nd:YAG  laser  receiv¬ 
er,  this  corresponds  to  a  relative  target  velocity  of  about 
5  mph  (2.1  mtrs/sec),  a  rather  small  velocity  and  doppler 
shift.  The  following  model  parameters  were  used  for  the 
simulation : 


I 

o 

o 

1 _ 

1 

o 

rH 

"3.6x10^^  1.09x10®' 

Q  = 

(t»  = 

Po  = 

_0  0.005. 

_T  1_ 

_1. 09x10®  ttVS  _ 

(95) 
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One  hundred  recursive  estimates  were  simulated  and  averaged 
over  20  realizations  of  the  measurement  noise. 

The  Po  matrix  was  chosen  iteratively  in  an  effort  to 
get  the  filter  to  converge  without  tracking  the  noise  (Ref. 
7:338).  Initially,  Poii  was  set  to  36  THz“  by  assuming 
that  the  frequency  error  or  doppler  is  normally  distributed 
with  standard  deviation  6  MHz.  The  phase  variance,  Poz2  , 
was  again  chosen  at  tt^/3  .  The  initial  trial  of  Poi2=Po2i 

was  /PoiiPo22=l. 09x10^  but  this  value  was  too  large  and 
caused  the  estimator  to  track  the  measurement  noise.  Ex¬ 
amination  of  the  estimator  equation  (57)  indicates  that  the 
element  influences  the  update  of  xi  ,  so  P012 

was  reduced  until  the  filter  converged  more  quickly  to  the 
actual  variance  of  its  estimates  (or  something  close  to  it 
in  this  case). 

At  30  and  lOdb  CNR,  the  filter  makes  the  following  av¬ 
erage  errors  and  standard  deviations: 

”.  56MHz  (020)11  =  .  357MHz 

CNR=30  db:  X2o~2i  “  (96) 

-.024rad_  (020)22  =  . 102rad 


CNR=10  db :  X20-X  = 


(020)11  =  .375MHz 

(97) 

(020)22  =  . 159rad 


The  processor  exhibits  little  degradation  in  performance 
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when  the  CNR  drops  from  30db  to  lOdb  compared  to  a  similar 
drop  when  the  frequency  is  known  exactly.  However,  the  var¬ 
iances  of  the  errors  are  significantly  higher  than  the  error 
variances  for  the  same  CNR's  when  the  frequency  is  known 
exactly.  No  comparisons  were  made  with  other  processors, 
but  the  phase  estimate  seems  better  than  might  intuitively 
be  expected  for  the  given  error  in  the  estimate  of  the  fre¬ 
quency  error  ("^30%).  This  is  not  that  surprising  since  the 
actual  frequency  is  52MHz,  an  error  in  the  estimate  of  the 
offset  of  500KI1Z  gives  a  fractional  frequency  error  of  less 
than  1%. 
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VI.  Conclusions  and  Recommendations 


A  nonlinear  Kalman . filter  has  been  derived  to  estimate 
the  phase  of  a  carrier  measured  through  an  IDF.  Two  models 
of  the  phase  process  were  examined  and  one  case  involving 
doppler  shift  was  also  examined,  though  not  in  depth.  The 
problem  was  motivated  by  an  optical  heterodyne  receiver 
structure,  but  the  solution  formulation  was  much  more  gen¬ 
eral  in  nature.  Any  system  which  processes  samples  like 
those  in  equation  (25)  could  use  all  or  part  of  the  theo¬ 
retical  development  in  this  thesis. 

Conclusions 

This  study  has  demonstrated  that  it  is  possible  to  es¬ 
timate  the  phase  of  an  optical  field  with  detectors  which 
are  much  lower  in  bandwidth  than  the  heterodyned  waveform 
at  the  IF.  The  only  requirement  placed  on  the  system  is 
that  the  integration  time  of  the  detector  be  small  enough 
that  the  phase  and  amplitude  can  be  treated  as  constants 
over  the  integration  period. 

Although  the  processor  shows  no  significant  improve¬ 
ment  over  the  linearized  PLL  for  the  two  phase  models  ex¬ 
amined,  the  PLL  could  be  utilized  only  in  conjunction  with 
wideband  detectors.  Thus,  this  system  allows  phase  measure¬ 
ments  where  they  would  otherwise  be  impossible.  Also,  for 


53 


cases  where  distributed  phase  measurements  are  desirable, 
wideband  detectors  are  generally  not  available  in  large  ar¬ 
rays,  but  CCD's  are. 

The  few  simulations  performed  have  also  indicated  that 
it  is  possible  to  estimate  relative  frequency,  doppler 
or  LO  instabilities,  in  addition  to  relative  phase.  These 
simulations  have  also  indicated  that  the  phase  can  be  track¬ 
ed  even  with  a  relatively  poor  estimate  of  the  frequency 
available.  The  system  is  not  too  sensitive  to  the  actual 
LO  frequency. 

Recommendations 

The  measurement  equation  (25)  was  derived  on  the  basis 
of  two  simplifying  assumptions:  the  amplitude  and  the  phase 
are  both  approximately  constant  over  the  measurement  inter¬ 
val.  Since  a  CCD  inherently  averages  over  a  rather  long 
period  of  time  compared  to  T  ,  the  first  assumption  would 
probably  be  violated  in  a  rapidly  fading  environment  or 
where  the  carrier  has  incidental  AM  of  wide  bandwidth  com¬ 
pared  to  the  integration  time.  Similarly,  the  second  as¬ 
sumption  could  be  investigated  by  allowing  the  phase  to 
drift  during  the  measurement  interva] . 

The  presence  of  a  modulo-ir  ambiguity  in  the  processor 
is  a  drawback,  but  it  might  be  possible  to  eliminate  it. 

One  possible  method  to  accomplish  this  was  suggested  by 
Meer  (Ref . 10: 106) .  That  method  consists  of  checking  to  see 
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if  0j^  is  close  to  ±  ^  .  If  it  is,  and  6)^+2 

close  to  ±  ^  and  has  sign  opposite  that  of  9^^ 
new  estimate  of  defined: 


is  also 
,  then  a 


(98) 


The  sign  in  front  of  the  tt  depends  on  whether  0j^  was 
positive  or  negative.  In  the  absence  of  noise,  this  pro¬ 
cedure  could  regenerate  the  phase  sequence  with  no  ambigu¬ 
ity.  In  the  presence  of  noise,  the  number  of  errors  com- 
mitted  would  depend  on  how  small  the  parameter 
is  permitted  to  become  and  how  small  I  ®k+l  ^  must 

alsc  "ecome  before  a  correction  is  made,  provided  9 
changes  sign.  Clearly,  more  errors  would  also  be  made  as 
the  CNR  is  decreased. 

This  procedure  is  somewhat  ad-hoc.  A  more  mathemati¬ 
cally  precise  method  would  be  to  take  advantage  of  both  the 
in  phase  and  quadrature  information  contained  in  the  measure¬ 
ment.  The  update  portion  of  the  covariance  matrix,  the  term 


ecos(  )+nj^^  in  equation  (66),  or  the  algebraically 

equivalent  part  of  equation  (58),  is  the  part  that  incorpor¬ 
ated  new  measurement  information  into  the  covariance  equa¬ 
tion.  If  l®k+l~®kl^^  ’  then  a  wrap  around  will  occur  and 

the  update  portion  of  the  variance  equation  will  be  negative, 
since  cos(x)<0  when  tt<|x|<2it  .  Adding  tt  to  9^+1 
it  falls  into  the  fourth  quadrant  (or  subtracting  tt  from 
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if  it  falls  into  the  first  quadrant)  when  the  update 
term  is  negative  will  regenerate  the  phase  sequence.  If  the 
noise  variance  becomes  too  large,  clearly  errors  may  be  made, 
but  the  estimator  does  not  perform  well  in  very  low  CNR  en¬ 
vironments,  anyway. 

Much  work  remains  open  on  the  problem  of  estimating  a 
multi-parameter  state  vector.  In  particular,  the  most  in¬ 
teresting  case  from  an  optical  communication  point  of  view 
is  the  unstable  LO  tracking  a  random  phase  walk.  The  fil¬ 
ter  equations  were  derived  in  (92)  and  (93)  (implicitly,  at 
least),  but  performance  simulations  were  only  begun  wdth  a 
small,  fixed  unknown  frequency  shift.  Extensions  of  this 
work  would  have  practical  applications  to  any  laboratory 
or  field  system  since  LO  instabilities  are  always  present. 

Since  the  main  advantage  of  CCD  detectors  is  their 
availability  in  large  arrays  which  allow  spatially  distri¬ 
buted  measurements,  it  would  be  logical  to  incorporate 
spatial  information  into  the  processor.  If  the  phase  of 
the  received  waveform  or  the  background  noise  has  a  known 
or  assumed  correlation,  then  the  entire  array  of  measure¬ 
ments  can  be  utilized  to  estimate  the  phase  at  each  pixel. 

The  incorporation  of  more  a-priori  information  should  re¬ 
sult  in  a  lower  variance  processor.  Because  of  the  conceiv¬ 
ably  large  size  of  the  array,  either  a  piecewise  correlation 
over  small  parts  of  the  array  would  be  assumed,  or  the  array 
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could  be  processed  in  a  spatially  recursive  fashion, 
second  alternative  is  especially  intriguing  since  the 
tion  of  the  permissibility  of  recursive  processing  in 
dimensions  is  not  quite  closed  (Ref , 11 ; 935 ) . 


The 

ques- 

two 
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Appendix  A 

Proof  that  Samples  of  a  Sine  Wave  Sum  to  Zero 


Numerous  times  in  the  text,  a  great  deal  of  mathematical 

simplification  resulted  when  terms  of  the  form  sin(^(i-l) 

Air  t 

+  0(t))  and  sin  i~l  )+o  ( t ) )  summed  from  i  =  l  to 

were  set  to  zero.  [lere  this  identity  is  established.  Con¬ 
sider  a  slightly  more  general  case  where  k  is  an  arbitrary 
integer : 


^  sin(|— k*  (i-l)  +  0(t))  =  ^  sin(2iri.p - +  0(t)) 

i=l  ^  i=l  ^  ^ 

(99) 

Euler's  identity  states  that  cxp( j x )=cos( x  )  + j s: n ( x  )  ,  so, 

if  it  can  bo  established  that: 


E 

i  =  l 


exp 


j{2.i 


,1^ 

'ic 


2^k 


0(1)} 


=  0 


(100) 


then  the  proof  is  completed.  The  term 
ca9  be  moved  out  of  the  summation  since 
on  the  index,  i  .  Define  a  parameter 
that  a  geometric  series  in  L  results, 
as : 


exp 


+jO(t) 


it  does  not  depend 
L=oxp( ,i2Trk/P^  )  so 


A  geometric  series 


sums 
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(101) 


t  p 

E  = 

i=l 

Pt 

Since  k  is  an  integer,  L  =0  unless  k=0  ,  so  the  sum 

in  equation  (100)  is  indeed  zero  for  all  k^O  ,  including 
tlie  special  cases  of  interest,  k=l  and  k=2 
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Appendix  B 

Simplification  of  Estimator  Equations 


The 

to 

For 


Several  vector  differentiation  rules  must  be  defined, 
operator  is  a  row  vector  which  is  always  multiplied 

the  riglit  of  tlie  vector  on  which  it  operates  (Ref.  1:983). 
a  3-dimensional  vector: 

6^(  x) 


f  1  (x) 

f  2  (x) 
f  3  (x) 


6  X  j  6x2 


&K-. 


(102) 


Similarly,  tlie  operator  is  a  column  vector  which  al¬ 

ways  operates  to  the  left  of  a  transposed  column  vector  (Ref. 


1:983) : 


6f^(x) 

6x 


6 

6x  1 


[fl(x)  f2(x)  f3(x)] 


(103) 


Both  differentiations  result  in  matrices  with  elements  like 

6f ^(x) 

6x  . 

J 

From  the  chain  rules  for  these  operators  (Ref. 1:985), 
the  partial  derivatives  in  equations  (49)  and  (50)  can  be  ex¬ 
panded  : 
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(104) 


5  ,T,  T,  -  , 

-k  6x,,  60' 


T  -^h^ep 

1>  £ - - 

60." 


(105) 


6  .T 


.T 

'^£k 


6h(ep  ^ 


6  0; 

k 


(106) 


Equation  (49)  becomes: 


_T  ^  ^®k^  ^-1 

60, 


Ip 


(107) 


<®p 


£k+l  "^-k  ^k+l£  ^  Ik+l"I^®k^] 

60k  L  J 


(108) 


which  is  equation  (53).  Equation  (50)  becomes: 


P  =  P' 
^k+1  k 


-T  -T 

<i>  -([> 

6  Tpc  «h’^(0k)  ij 

Ik+l-I(®P 

-1 

T  ^ 

r-  '  P 

£  ^k 

60j^L  60j^ 

-1 

- 1 


(109) 


=  p; 


, -T 


f 

"«8;  L 


^  /N 

6h^0p  _, 

- —  R 


,  i_  so; 

k  I-  k 


Ik+l-I^Op 


T 

£  Pk 


- 1 


(110) 
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The  vector  operators  also  obey  the  rule  (Ref. 1:985): 

r  A  A  *5  Z 

_(Az)  =  _  .  z  +  A  .  ^  (111) 

T 

5h^(en 

let  A= -  and  let  z=R  '{r.  ^-h(9')}  .  Then: 

69"  -K+i  K 


P  =  P' 
k+1  k 


-T 

«"htep 

y\ 

$  -c 

- p 

^  'a  ^ 

60- 

r,  ,  1  -h( 0,''  ) 
—k+1  —  k 

Sh'^COn  6h(0p 

- *i_  .  R  - JL- 


T  ^ 

^  Pk 


_  1 


(112) 


Which  is  equation  (54). 

Now  substitute  equations  (25)  through  (27)  into  equa¬ 
tions  (53)  and  (54)  along  with  R=a^I  ,  where  I  is  a 

n 

P^xP^  identity  matrix.  Write  out  the  vector  products  as 
summations  to  give; 


-  -  1  ^ 

X,  ,  -  =  (px,  +  P,  c  •  -  - 

—k+1  k  k+1  —  Z_/  IT 


sin(p:^)sin(p^(i-l  )  +  0p 


^k+l^i) 


i  Ty 

;  sin(^)cos(^(i-l)+0j;) 


i  T 
s 


(113) 


Utilizing  Appendix  A 

-  - 

X,  .  1  =  <px,  -  - z-’sin 

—k+1  k  TTO 


and  some  trigonometric  identities; 


(114) 
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which  is  equation  (57).  Before  doing  similarly  for  equa¬ 


tion  (54),  move  the  inside  the  inversion  operation  on 

the  brackets: 


k+1 


-T  ^-1 
i  Pk  -£ 


6^h'(0p  , 

- —  R~' 

60- 


60, 


6  0; 


- 1 


(115) 


where  the  fact  that  P^Pr”'^!  utilized.  Now  let  R  ^  = 

l/o^I  : 

'  n 


k+1 


-Tr,.-i  1 


6’h'‘'<6;) 


%(0k> 

2 

T 

60^ 

c 

- 1 


(116) 


Substituting  for  P'  from  equation  (55)  and  utilizing  (25)- 
(27): 

P. 


k+1 


T  T  - 1  T  1 

((J,Pj^(^'+GQG')  -  c  c"  ^ 


n 


i  Ty 

-f-  • 


•  cos(^(i-l)+0p  *'k+l(^)" 


i  Ty 

sin(~)cos(|^(i-l) 
t  t 


i  Tt  i  Ty 

_ § _ V(-^ _ 

P^.  77 

Z  J 


^  )  sin^  (^)sin^(|-^(i-l)  +  0pi 


- 1 


(117) 
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The  defined  in  equation  (63)  appears  naturally. 

Now  substitute  this  into  equation  (114),  use  Appendix  A, 

i  Y 

and  multiply  by  (-^)(j^)  =  l  . 

s  ' 


i  TY 

j,  s 

= 


i  Y 


— k+1  k  TTO 


n 


t  s 


■sin(^(i-l)  +  0j^) 


(124) 


1 


.  /S  ,  X  ..TT.  TT  V— V  ,  .  . 

(px,  -  ( - )  • — r*sin(K— )*-^ -  >  z,  ,,(i) 

k  TT  '  ^t  ^s^  ^  k+1^ 


Ptt  a 

.sin(^(i-l)  +  6p 


(125) 


=  4iXi,  - 


s 


sin(^) 


(126) 


2P. 


=  - 


-k  i  Ytt 
s 


CNR 


•sin(^)2^Zj^^^(i)sin(|^(i-l)  +  0') 


(127) 


The  factor  CNR  was  defined  in  (62)  and  can  be  verified  by 
substituting  equation  (15)  for  and  remembering  that: 


n'(i)  -  i  n^(l)  *  (,e  .  ^  a‘ 


(128) 
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appears  instead  of 


In  equation  (60)  the  factor  —  appears  instead  of 

s’ 

because  it  is  assumed  throughout  the  paper  that  the  signal 
power  is  normalized,  4(1^^)^=!  ,  and  fluctuations  in  the 


CNR  are  modelled  only  as  fluctuations  in  o 


,  0  ,  or 

n  ’  n'  ’ 


No  .  The  exact  same  arguments  for  equation  (58)  would  lead 
to  equation  (61). 

Now  examine  the  term  in  equation  (127); 


-|-sin(^)cos(~(i-l)  +  ej,^^) 


'sin(|^(i-l)  +  0p+n'^^(i)sin(|l-(i-l)+0p 


(129) 


=  ^“§7  sin(~)sin(0'-0j^^^)+n^^j^(i)sin(|^(i-l)  +  Op(13O) 


1  TP. 

-  sin(^)  sin(0j,^^-0p  +  nj^^ 

t 


(131) 


Similarly  for  equation  (118): 


(132) 
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whore  n,  is  defined  in  an  obvious  manner.  Since  n, 

ks 

and  nj^^  are  both  the  sums  of  weighted  gaussian  variables, 
they  are  also  gaussian,  so  it  suffices  to  compute  their 


second  moments. 


[■’ks] '  *=  E" 


'^^(i)sin(|^(i-l)+0p 


(133) 


E^e  E„["ktiO)]sin(|l(i-i)k5') 


(134) 


(135) 


(136) 


(137) 


Where  the  identity  j^x(y)j  =Ey  h|yH^]]  was  used,  as 

was  the  fact  that  *^^(1)  Is  zero  mean.  It  is  easily  shown 
that  ^Kc]=''  •  s  ince  nj^^  and  nj^^  are  zero  mean, 

their  variances  are  their  second  moments. 


ks"k's]  =  ® 

_  1  J 


(k+l)  n'(k+l)sin(|^(i-l)+e') 


Ptt  ^ 

•sin(p-(  j-l)  +  0j^. 


(138) 
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= 


EE  E  [nr(k+l)n-:(k+l)]- {cos(|^-(i-l) 

ijOLi  j  j 


/V 


.2-n 


iA,  ^ 


+  0j^-e'.  )“COs(p-(  i+j-2)  +  9j^+0j^.  ) } 


(139) 


=  E; 


2  E^[n^(k+l)n^(k4:i)]  -i 


(140) 


kk' 


(141) 


The  fact  that  j^n((k+l)nj(k'+l)j  was  used  since 

n/”  is  white.  Similarly, 


E 


E 


=  0 


(142) 


(343) 


So,  spectrally  white  seqi’  ..os,  at- 

ually  independent,  and  gaussian. 

The  above  derivations,  plus  letting  (iY)^/2=l  ,  re- 

O 

suit  in  equations  (65)  and  (66). 
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Appendix  C 

Simulation  Subroutine 
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n  n  n 


SUBRDUT I NE  GENS  I M ( Y , P , YO , PO , YACT , Q , PH  I , K  MAX , I  TER . 
•tlYDIM,  I  PT.CNRDB,  SUMP,  SUMY,  N,KMPT,WK  AREA) 

DOUBLE  DSEED 
REAL  N 

DIMENSION  Y( lYDIM.KMAX) ,P( lYDIM, lYDIM.KMAX) , YO( lYDIM) 
DIMENSION  PCX  lYDIM, lYDIM) ,Q( lYDIM, lYDIM) 

DIMENSION  YACT(IYDIM,KMAX) , PH I ( I  YD  I M , I  YD  I M ) 

DIMENSION  SUMY( IYDIM,KMAX) ,SUMP( lYDIM, lYDIMrKMAX) 
DIMENSION  WKAREACKMPT) ,N(KMPT) 

C 

c 

C  2D  LT  MARTIN  B.  MARK,  GEO-BOD,  22  SEP  80 
C 

C  THIS  SUBROUTINE  IS  A  GENERALIZED  SIMULATOTy  FOR  THE  DE- 
C  MODULATION  OF  A  SIGNAL  ENCODED  ON  THE  PHASE  OF  A  RINE- 
C  WAUE.  VIRTUALLY  ANY  PHASE  PLANT  MODEL  CAN  BE  ACCOMD- 
C  DATED.  INTEGRATED  SAMPLES  OF  TflE  RECEIVED  WAVEFORM 
C  ARE  GENERATED  AND  USED  IN  AN  EXTENDED  KALMAN  FILTER 
C  TO  DEMODULATE  THE  SIGNAL. 

C  THE  VARIABLES  USED  HAVE  THE  FOLLOWING  SIGNIFICANCES.’ 

C 

C  Y( I ,K  ) 

C  P<I,J,K) 

C  YO(I) 

C  PO<I,J) 

C  YACT<I,K) 

C  Q  (  I  ,  J  ) 

C  PHI(I,J) 

C  KMAX 

C  ITER 

C 

C  lYDIM 

C  IPT 

C 

C  CNRDB 

C  KMPT 

C 

C  SUMP, SUMY, 

C  N,WKAREA 

C 

C  IMSL  LIBRARY  ROUTINES  GGNML  AND  LINV3F  ARE  USED  TO 
C  GENERATE  GAUSSIAN  NOISE  SAMPLES  AND  INVERT  A  MATRIX, 

C  RESPECTIVELY.  SUBROUTINES  NEWP,  NEWY,  AND  SUMM  ARE 
C  ALSO  UTILIZED  AND  ARE  ATTACHED  AT  THE  END  OF  GENSIM. 

C  THE  PHASE  ANGLE  IS  ALWAYS  THE  LAST  ELEMENT  OF  THE 
C  STATE  VECTOR  Y.  ALL  OTHER  ELEMENTS  OF  THE  STATE  VEC- 
C  TOR  ARE  ALSO  COMPUTED  AND  RETURNED,  BUT  THEIR  PHY- 
C  SICAL  SIGNIFICANCES  ARE  ARBITRARY. 

C 

CHUHHf*#*****#****##****#*  **##*#■»******■»«■■»»«••»  ****«•«*«»**»•»■»# 
INITIALIZE  CONSTANTS 


K-TH  STATE  VECTOR  ESTIMATE 

K-TH  VARIANCE  MATRIX 

INITIAL  STATE  VECTOR 

INITIAL  VARIANCE  MATRIX 

ACTUAL  STATE  VECTOR  SEQUENCE 

PLANT  NOISE  VARIANCE  MATRIX 

PLANT  TRANSITION  MATRIX 

NUMBER  OF  SEQUENTIAL  ESTIMATES  OF  Y 

NUMBER  OF  ITERATIONS  Y  IS  AVERAGED 

OVER 

DIMENSIONALITY  OF  STATE  VECTOR 
SAMPLES  OF  INPUT  WAVEFORM  PER 
PERIOD 

CARRIER-TO-NOISE  RATIO,  DB 
NUMBER  OF  NOISE  SAMPLES  TO  BE  GEN¬ 
ERATED  BY  GGNML 

SCRATCHPAD  MATRICES 
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KWD=2*IYDIK 
PI=3. 141592G54 
PT=FLOAT( IPT) 

CNR=10.**(CNRDB/10. ) 

ALPHA  =  SQRT(2. ) *CNR»PT*S I N ( P I / PT ) / P I 
D9EED=1073245SBO.DOO 

CLEAR  SUMP  AND  SUMY  MATRICES, 

SET  FIRST  ESTIMATE  OF  Y  TO  YO  AND 
FIRST  VARIANCE  MATRIX  TO  PO 

DO  120  I=1,IYDIM 
DO  no  J  =  1,IYDIM 
DO  100  K  =  1  rKMAX 
SUMP< I . J . K ) =0.0 
SUNY( J . K ) =0.0 
100  CONTINUE 

P( I , J, 1 )=P0( I, J) 
no  CONTINUE 
Y( I , 1 ) =Y0(  I  ) 

120  CONTINUE 

AVERAGE  QUER  "ITER"  REALIZATIONS  OF  THE  NOISE 
PROCESS  BY  INCREMENTING  DSEED  ON  EACH  PASS. 
NOISE  SAMPLES  ARE  GENERATED  BY  SUBROUTINE 
GGNML. 

DO  500  1=1, ITER 

CALL  GGNML(DSEED,KMPT,N) 

COMPUTE  SEQUENCES  OF  Y  AND  P  FOR  ALL  K 


KMAX=KMAX-1 
DO  400  K=1 rKMAX 
KK  =  IPT*<K-1  ) 

COMPUTE  A-PRIORI  UPDATE  OF  STATE  VECTOR 

ARGSUP=0.0 

DO  210  IND=1 , lYDIM 

ARGSUP=PHI ( I YDIM, IND)*Y ( INDrK ) +ARGSUP 
210  CONTINUE 

SUMS  AND  SUMC  ARE  THE  CORRELATIONS  BETWEEN  Z 
AND  THE  SINE  AND  COSINE,  RESPECTIVELY,  OF  THE 
A-PRIORI  ESTIMATE  OF  THE  PHASE. 


BUMC=0.0 
SUMS=0.0 
DO  300  J=1 , IPT 

N(KK  +  J  )  =N(KK-*-J  )/SORT(CNRi)PT> 

CON  =  2  .  ^^PI  «FLCAT  (  J-1  )  /  PT 
ARG=CON+ARGEUP 

Z= < ALPHA/ <CNR*PT) ) «COS ( CON+ YACT ( IYDIM,K) )+ 
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*  N(KK+J) 

SUMS  =  Z<^SIN(ARG)+S'JMS 
SUMC=2*C05 ( ARG)+SUMC 
300  CONTINUE 

C 

C  UPDATE  Y  AND  P  USING  SIMULATED  DATA  AND  PLANT 
MODEL. 

CALL  NEWP<  P, PO, PHI ,Q, I  YD  I M , SUMC , ALPHA , K , KMAX , 
«  WKAREA,KND,N) 

CALL  NEWY( Y, PHI .P. I  YD I M, SUMS , ALPHA , K ,KMAX, PI ) 
CALL  SUMM(£UMY,SUMP,YrP, I YDIM , K , KMAX ) 

400  CONTINUE 

INCREMENT  DSEED. 


KMAX=KMAX+1 
DSEED =DSEED+1 00. DOO 
500  CONTINUE 

COMPUTE  AOERAGE  Y  AND  P  FROM  SUMMATION  OF  ALL 
"ITER"  REALIZATION  OF  NOISE  PROCESS. 

DO  G20  K=1 rKMAX 

DO  610  J=1 , lYDIM 
DO  GOO  1  =  1 , I  YD  I M 

P(I , J ,K ) =SUMP( I , J,K) /FLOAT! ITER) 

GOO  CONTINUE 

Y ( J , K ) =SUMY  <  J , K ) /FLOAT ( ITER ) 

GIO  CONTINUE 
ezo  CONTINUE 


RETURN 

END 
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SUBROUTINE  NENP ( P , PO , PHI ,Q, I  YD  I M , SUMO , ALPHA . K , KMAX , 
■#'NKARtAfKND»S) 

DIMENSION  P( lYDIM, lYDIM.KMAX) , PHI ( I YD  I M , I YD I M ) 
DIMENSION  0( lYDIM, lYDlM) , PO < I  YD I M , I  YD  I M ) 

DIMENSION  WKAREA(KND) 

C 

C  SUBROUTINE  NEWP  COMPUTES  THE  NEXT  P  GIUEN  THE 
C  LAST  P,  SUMS,  AND  THE  PLANT  MODEL.  SUBROUTINE 
C  LINU3F  IS  USED  TO  INUERT  MATRICES  WHEN  lYDIM 
C  IS  NOT  EQUAL  TO  ONE. 

C 

DO  400  1=1 , lYDIM 
DO  300  J=1 , lYDIM 
DUNMYP=0.0 

DO  200  M=1.IYDIM 

DO  100  N=1 , lYDIM 

DUMMYP  =  PHI  <  I  ,N>*P(N,M,K  )#P1II  (  J  ,M)+DUMMYP 
100  CONTINUE 

200  CONTINUE 

PO(I , J ) =DUMMYP  +  Q( I , J  ) 

300  CONTINUE 
400  CONTINUE 

IF( lYDIM.EO. 1 )  GOTO  720 
Dl  =  l  . 

CALL  LIN03F(P0,S- 1 , I  YD IM , I YD I M , D 1 , D2 , WK AREA , I ER ) 

IF( lER.EQ. 130)GGTQ  700 

500  POdYDIM,  lYDIM)  =  P0<  lYDIM,  IYDIM)+ALPHA*SUMC 

CALL  LIN'.'3F(  PO,S,  1  ,  I  YD  I M  d  YD  I M  ,  D 1  ,D2,WKAREA,  lER) 

IF( lER.EQ. 130)G0T0  710 
510  DO  GIO  I  =  j  ,  lYDIM 

DO  GOO  J= 1 , 1 YDIM 
P( I  f  J , K  +  1 ) =PC ( I r  J  ) 

COO  CONTINUE 

GIO  CONTINUE 


700 


C 

710 


C 

C 

rj 


720 


RETURN 

PRINT  "SINGULAR  MATRIX  ENCOUTERED  IN  NEWP" 
PRINT  _LED  BY  FIRST  OCCURENCE  OF  LIN03F" 

GOTO  500 

PRINT  "SINGULAR  MATRIX  ENCOUNTERED  IN  NEWP" 
PRINT  *,  "CALLED  BY  SECOND  OCCURENCE  OF  LIN03F" 
GOTO  ;>io 


P( 1 , 1 ,K+1 ) = 1 . / ( ( 1 . /P0( 1 d ) >+ALPHA#SUMC) 


RETURN 

END 
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SUBROUTINE  NEWY ( Y , PH  I , P , I  YD  I M , SUMS , ALPHA , K , KMAX , P I ) 
DIMENSION  Y( lYDIM, KMAK) . PH I ( I  YD  I M , I  YD  I M ) 

DIMENSION  P{ lYDIM. IYDIM,KMAX) 

C 

C  SUBROUTINE  NEUY  COMPUTES  THE  NEXT  Y  GIUEN  THE 
C  LAST  Y,  THE  NEW  P,  SUMS,  AND  THE  PLANT  MODEL.  THE 
C  PHASE,  Y{ IYDIM,K+1 ) ,  IS  RETURNED  M0DUL0-2«PI 
C 

DO  200  I=1,IYDIM 
DUMMYY=0.0 

DO  100  J=1,IYDIM 
DUMMYY=PHI ( I , J )#Y( J,K )+DUMMYY 
100  CONTINUE 

Y( 1 ,K+1 ) =DUMMYY-P< I , IYDIM,K+1 >«ALPHA*SUMS 
200  CONTINUE 

Y(  IYDIM,K+1 ) =AMOD(Y( IYDIM.K+1 ) , PI ) 


RETURN 

END 


SUBROUT  I NE  SUMM ( SUMY , SUMP , Y , P , I  YD  I M , K , K  MAX ) 
DIMENSION  SUMY( IYDIM,KMAX) ,SUMP( lYDIM, IYDIM,KMAX) 
DIMENSION  Y( IYDIM,KMAX) , P( lYDIM, lYDIMrKMAX) 

SUBROUTINE  SUMM  ACCUMULATES  A  SUMMATION  OF  P 
AND  Y  FUR  ALL  "ITER"  REALIZATIONS  OF  THE  NOISE 
PROCESS. 

100  DO  300  1=1 , lYDIM 

DO  200  J= 1 , lYDIM 

SUMP( I , J ,K ) =P< I , J ,K)+SUMP( I , J , K ) 

200  CONTINUE 

SUMY( I , K ) =Y< I ,K )+SUMY( I ,K ) 

300  CONTINUE 

IF(K .LT.KMAX)  RETURN 
IFU.GT.KMAX)  GOTO  400 
K  =  K  +  1 
GOTO  100 

400  K=KMAX 


RETURN 

END 
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